/*
 * fft.cpp
 *
 *  Created on: 06.02.2013
 *      Author: Roma Jam
 */
#include "fft.h"

using namespace std;
fft_const::fft_const (int inv){
    int l, le;

    for (l = 0, le = 1; l < 15; l++, le <<= 1){
        if (inv) c[l] = exp((complex<double>(0.0,1.0)) * M_PI/(double)le);
        else     c[l] = exp((complex<double>(0.0,-1.0)) * M_PI/(double)le);
    }
} // fft_const

static fft_const fftcnst(0);
static fft_const ifftcnst(1);

complex<double>* fft (complex<double>*d, complex<double> *sa, complex<double> *sb, int n, int inv){
    complex<double> u,t;
    int nv2,nm1,i,j,k,l,le,le1,ip,m;

    for(m = 0,l = 1; l < n; m++, l <<= 1);  // m = log2(n)
    nv2 = n >> 1;
    nm1 = n-1;

//    if ( d != s ){
        for (i = 0; i < nv2; i++) d[i] = sa[i];
        for (i = nv2,j=0; i < n; i++,j++) d[i] = sb[j];
//    }
    for (j = 0, i = 0; i < nm1; i++, j += k){
        if (i < j){
            t = d[j];
            d[j] = d[i];
            d[i] = t;
        }
        for ( k = nv2; k <= j; j -= k, k >>= 1 );
    }
    for (l = 0, le1 = 1; l < m; l++, le1 <<= 1){
        le = le1 << 1;
        u = complex<double>(1,0);

        for ( j = 0; j < le1; j++ ){
            i = j;
            do{
                ip = i + le1;
                t = u * d[ip];
                d[ip] = d[i] - t;
                d[i] += t;
                if(inv){
                    d[ip] /= 2.0;
                    d[i] /= 2.0;
                }
                i += le;
            } while ( i < n );
            if(inv) u *= ifftcnst.c[l];
            else    u *= fftcnst.c[l];
        }
    }
  return d;
}

complex<double>* fftm (complex<double> *d, int n, int inv){
    complex<double> u,t;
    int nv2,nm1,i,j,k,l,le,le1,ip,m;

    for( m = 0, l = 1; l < n; m++, l <<= 1 );   // m = log2(n)
    nv2 = n >> 1;                               // nv2 = n / 2;
    nm1 = n - 1;                                // nm1 = n - 1;

//    if ( d != s )
//    for ( i = 0; i < n; i++ ) d[i] = s[i];

    for(j = 0, i = 0; i < nm1; i++, j += k){
        if ( i < j ){
            t = d[j];
            d[j] = d[i];
            d[i] = t;
        }
        for(k = nv2; k <= j; j -= k, k >>= 1);
    }
    for(l = 0, le1 = 1; l < m; l++, le1 <<= 1){
        le = le1 << 1;
        u = complex<double>(1,0);
        for(j = 0; j < le1; j++){
            i = j;
            do {
                ip = i + le1;
                t = d[ip] * u;
                d[ip] = d[i] - t;
                d[i] += t;
                if (inv){
                    d[ip] /= 2.0;
                    d[i] /= 2.0;
                };
                i += le;
            } while ( i < n );
            if(inv) u *= ifftcnst.c[l];
            else    u *= fftcnst.c[l];

        }
    }
  return d;
}
